---
title: "CR_Part1b_compileMLmodANOVAtab"
author: "brouwern@gmail.com"
date: "January 9, 2017"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


# Intro

Preceding file as of 1/9/2017 is:

"Latta_et_al_2017_PeerJ_multilevel_model_with_predictors.Rmd"

This file runs the multilevel pois-normal mixed model with different traits as predictors

This file processes the output of that, specifically it

* extract ANOVA tables that were compiled prevviously and put into single table
* prep table for output to Excel; most labeling is done here


* Compile ANOVA table that becomes Supp Table 1 (?)
* Make table of regression models; becomes Supp Table 3 (?)
* Plot MARGINAL effects FIgure 3 (?)
* Calculate "real"" slope for each group: Supp Table 2 (?)





# Re-load data
# Re-load all model output

All data was saved in lists of processed models and processed output
```{r}
#file names of output from modeling file
all.output <- c("canopy_ob_output_list.RData"
,"canopy_obYN_list.RData"

,"canopy_use_output_list.RData"
,"canopy_useYN_list.RData"

,"cons_prior_model_list.RData"
,"cons_priority_output_list.RData"

,"dist_sense_model_list.RData"
,"dist_sense_output_list.RData"

,"elev_mig_model_list.RData"
,"elev_mig_output_list.RData"

,"forage_guilde2lev_list.RData"
,"forage_guilde2lev_output_list.RData"

,"habitat_breadt_output_list.RData"
,"hab.breadth_list.RData"

,"forage_guilde3lev_list.RData"
,"forage_guilde3lev_output_list.RData"

,"hab_prime_model_list.RData"
,"hab_prime_output_list.RData")


#load all output into memory
for(i in 1:length(all.output)){
  file.i <- paste("./MODEL_OUTPUT/",all.output[i],sep = "")
  load(file.i)
  
}



```




# Compile ANOVA tables

rbind() all of the seperate ANOVA tables
```{r}

#for spacing between each table
blank.x <- rep(NA,dim(hab.prim.output.list$anova.table)[2])

#compile each table w/rbind()
#NB: each item is a serpate previously compiled table
#this just stacks all of the tables
all.anova.tables <- rbind(hab.prim.output.list$anova.table,blank.x
                          ,dist.sense.output.list$anova.table,blank.x
                          ,cons.priority.output.list$anova.table,blank.x
                          ,elev.mig.output.list$anova.table,blank.x
                          ,canopy.ob.YN.output.list$anova.table,blank.x
                          ,canopy.use.YN.output.list$anova.table,blank.x
                          ,forage.guilde.2lev.output.list$anova.table, blank.x
                          ,forage.guilde.3lev.output.list$anova.table, blank.x,
                          hab.breadth.output.list$anova.table)
```


clean things up 
```{r}
#names of each variable
model.names <- c("hab.prime"
                 ,"sense2"
                 ,"cons.priority"
                 ,"mig.elevYN"
                 ,"canopy.ob.YN"
                 ,"canopy.use.YN"
                 ,"forage.guilde.2lev"
                 ,"forage.guilde2"
                 ,"hab.breadth")

d <- dim(hab.prim.output.list$anova.table)[1]
repnames <- function(x,d = 7){rep(x,d)}

names.column <- NA
for(i in 1:length(model.names)){
  names.column <- c(names.column,
                    repnames(model.names[i], d = d),NA)
}
names.column <- names.column[-1]
names.column <- names.column[-length(names.column)]



all.anova.tables2 <- data.frame(trait = names.column,
           all.anova.tables)

names(all.anova.tables2)[2] <- "Test"

all.anova.tables2$df <- with(all.anova.tables2,
                             paste(df1,df2, sep = ","))

all.anova.tables2$"" <- ifelse(all.anova.tables2$p > 0.11,"",".")
all.anova.tables2[,8] <- ifelse(all.anova.tables2$p > 0.051,
                               all.anova.tables2[,8],
                               "*")

all.anova.tables2$trait <- as.character(all.anova.tables2$trait)

```



Make some labels for each sep. ANOVA table
```{r}


seq.i <- seq(1,dim(all.anova.tables2)[1], by = 8)
all.anova.tables2$trait[seq.i[1]] <- "Habitat preference (secondary vs. primary forest)"
all.anova.tables2$trait[seq.i[2]] <- "Sensitivty to disturbance (low vs. med/high)"
all.anova.tables2$trait[seq.i[3]] <- "Conservation priority" 
all.anova.tables2$trait[seq.i[4]] <- "Elevational migrant"
all.anova.tables2$trait[seq.i[5]] <- "Canopy Use (Facultative vs. Obligate)"
all.anova.tables2$trait[seq.i[6]] <- "Canopy Use (Use vs. Non-use)"
all.anova.tables2$trait[seq.i[7]] <- "Foraging guild (2 Levels: Omnivore vs. Specialist)"
all.anova.tables2$trait[seq.i[8]] <- "Foraging guild (3 Levels: Omnivore, Frug/Nectarivore, Insectivore)"
all.anova.tables2$trait[seq.i[9]] <- "Habitat breadth (Continous variable)"



all.anova.tables.fin <- all.anova.tables2[,c(1,2,5,7,6,8)]
round_df(all.anova.tables.fin,4)



```

OUtput the table as csv
```{r}
write.csv(all.anova.tables.fin, file = "ANOVA_table_all_multilevel_models2.csv")

```



# Make table of full regression models 

This is Supplemental table 1 and shows original, unprocessed regression coefficients.                 
                 
                 
## Function to extract coefficients from model list                 
```{r}
get.betas <- function(mod.list){
  betas <- summary(mod.list$m7a.best)$coefficients[,"Estimate"]
  ses <- summary(mod.list$m7a.best)$coefficients[,"Std. Error"]
  
  #ses <- paste("(",round(ses,2),")",sep = "")
  
  out <- list(betas = betas,SEs =ses)
  
  
  return(out)
}



```


Extract betas and SEs from a model from the list
```{r}
hab.prime.beta <- get.betas(hab.prime.mod.list)
dist.sense.beta <- get.betas(dist.sense.mod.list) #sense2
cons.prior.beta <- get.betas(cons.priority.mod.list)#cons.priority
elev.mig.beta <- get.betas(elev.mig.mod.list) #mig.elevYN
canopy.ob.YN.beta <- get.betas(canopy.ob.YN.mod.list) #canopy.ob.YN
canopy.use.YN.beta <- get.betas(canopy.use.YN.mod.list)#canopy.use.YN
forag.guilde.2lev.beta <- get.betas(forage.guilde.2lev.mod.list)  #forage.guilde.2lev
forage.guilde.3lev.beta <- get.betas(forage.guilde.3lev.mod.list) #forage.guilde2
hab.breadth.beta <- get.betas(hab.breadth.mod.list) #hab.breadth


```



Set up betas for table
```{r}
betas <- rbind(hab.prime.beta$betas
,dist.sense.beta$betas
,cons.prior.beta$betas
,elev.mig.beta$betas
,canopy.ob.YN.beta$betas
,canopy.use.YN.beta$betas
,forag.guilde.2lev.beta$betas
#,forage.guilde.3lev.beta$betas
,hab.breadth.beta$betas)


#names of each variable
model.names <- c("hab.prime"
                 ,"sense2"
                 ,"cons.priority"
                 ,"mig.elevYN"
                 ,"canopy.ob.YN"
                 ,"canopy.use.YN"
                 ,"forage.guilde.2lev"
                 ,"forage.guilde2"
                 ,"hab.breadth")

#drop forage.guilde2 b/c it has 3 factor levels
betas2 <- data.frame("Trait Name" = model.names[-8],betas)

names(betas2)[2] <- "Intercept"
names(betas2)[3] <- "Beta.trait"
names(betas2)[4] <- "Beta.Month.January"
names(betas2)[5] <- "Beta*Year"
names(betas2)[6] <- "Beta*Year*Trait"

```




Set up SEs

```{r}
SEs <- rbind(hab.prime.beta$SEs
,dist.sense.beta$SEs
,cons.prior.beta$SEs
,elev.mig.beta$SEs
,canopy.ob.YN.beta$SEs
,canopy.use.YN.beta$SEs
,forag.guilde.2lev.beta$SEs
#,forage.guilde.3lev.beta$SEs
,hab.breadth.beta$SEs)


SEs2 <- data.frame(SEs)

names(SEs2)[1] <- "a"
names(SEs2)[2] <- "b"
names(SEs2)[3] <- "c"
names(SEs2)[4] <- "d"
names(SEs2)[5] <- "e"



```



Combine betas and SEs 
```{r}

coefs.with.ses <- data.frame(Trait.name = betas2[,1]
,Intercept = betas2[,2],
SE = SEs2[,1]
,Beta.trait = betas2[,3],
SE = SEs2[,2]
,Beta.Month.January = betas2[,4],
SE = SEs2[,3]
,"Beta*Year" = betas2[,5],
SE = SEs2[,4]
,"Beta*Year*Trait" = betas2[,6],
SE = SEs2[,5])
```


Format the table

```{r}
names(coefs.with.ses)[1] <- "Trait.name"
coefs.with.ses$Model <- c(
  "Habitat preference",
  "Sensitivity to disturbance",
  "Conservation priority",
  "Elevational migrant",
  "Canopy Use-Obligate?",
  "Canopy Use-Any?",
  "Foraging  guild-2",
  #"Foraging  guild-3",
  "Habitat breadth ")

```



Save the table as csv
```{r}
write.csv(coefs.with.ses,file = "./model_output/regression_table.csv")
```





# Plot MARGINAL effects

These are plots of MARGINAL effects (NOT what I want in final paper!)

```{r}
ggplot(data = coefs.with.ses,
       aes(y = exp(Beta.Year.Trait),
           x = Trait.name)) +
  geom_point(size = 5)  +
  geom_errorbar(aes(ymax = exp(Beta.Year.Trait+ 1.96*SE.4),
                    ymin = exp(Beta.Year.Trait- 1.96* SE.4)),
                width = 0, size = 1.2) +
  theme_bw() +coord_flip() +
  geom_hline(yintercept = 1)
```






# Calculate "real"" slope for each group: Supp Table 2

* THis is what end up in Supp Table 2
* This is the "real" slope of the trend, not the marginal effect/beta

## Use multcomp to combine parameters of model

```{r}


calc.realized.slopes <- function(working.list){
library(multcomp)
working.model <- working.list$all.models$m7a.best
contrast.mat <-  t(as.matrix(c(0, 0, 0, 1, 1)))

contrast.mat <-  rbind(c(0, 0, 0, 1, 0),
                       c(0, 0, 0, 1, 1))


levs <- levels(working.model@frame$trait)

rownames(contrast.mat) <- levs

summary(working.model)$coefficients

glht.out <- glht(working.model,
     linfct = contrast.mat,
     alternative = c("two.sided"))

glht.out2 <- summary(glht.out,test = adjusted(type = "none"))



trait.names <- paste(working.list$trait,levs,sep = ".")
df.out <- data.frame(Trait = trait.names,
           Effect.size = glht.out2$test$coefficients,
           SE = glht.out2$test$sigma,
           p = glht.out2$test$pvalues)

df.out$sig <- ifelse(df.out$p > 0.2, "","-")
df.out$sig <- ifelse(df.out$p < 0.10, ".",df.out$sig)
df.out$sig <- ifelse(df.out$p < 0.05, "*",df.out$sig)
return(df.out)
 
}
  



```



  

### Make a table
Calcualte realized slopes and put into a table 

Conservation priority: 1 = urgent; 2 = high; 3 = medium; 4 = low


Extract and calc slopes for each sep group
```{r}
hp <- calc.realized.slopes(hab.prim.output.list)
ds <- calc.realized.slopes(dist.sense.output.list)
cp <- calc.realized.slopes(cons.priority.output.list)
em <- calc.realized.slopes(elev.mig.output.list)
c1 <- calc.realized.slopes(canopy.ob.YN.output.list)
c2 <- calc.realized.slopes(canopy.use.YN.output.list)
fg <- calc.realized.slopes(forage.guilde.2lev.output.list)
#calc.realized.slopes(forage.guilde.3lev.output.list)
#calc.realized.slopes(hab.breadth.output.list)





```


Label the tabel
```{r}
hp$Trait <- "Habitat preference\n(Primary vs.\nSecondary Forest)"
hp$Group1 <- c("Primary forest","Secondary forest")
hp$Group <- c("Mature forest","Early successional")

ds$Trait <- "Disturbance sensitivity\n(High/Med. vs.\nLow sensitivity)"
ds$Group1 <- c("High/Med","Low")
ds$Group <- c("Mature forest","Early successional")

cp$Trait <- "Conservation priority\n(Medium vs.\nLow priority)"
cp$Group1 <- c("Medium","Low")
cp$Group <- c("Mature forest","Early successional")

em$Trait <- "Elevational migrant\n(Migrant vs.\nNot elev. mig.)"
em$Group1 <- c("Elev. mig","Not elev. mig")
em$Group <- c("Mature forest","Early successional")

c1$Trait <- "Canopy use\n(Obligate vs.\nFac./Not used)"
c1$Group1 <- c("Obligate","Facultative/Not used")
c1$Group <- c("Mature forest","Early successional")

c2$Trait <- "Canopy use\n(Obligate/Fac. vs.\nNot used)"
c2$Group1 <- c("Oblig./Fac.","Not used")
c2$Group <- c("Mature forest","Early successional")

fg$Trait <- "Foraging guild\n(Omnivore vs.\nSpecialist)"
fg$Group1 <- c("Omnivore","Specialist")
fg <- fg[c(2,1),]
fg$Group <- c("Mature forest","Early successional")



```


Load object for latitudinal migration

```{r}
load(file = "./MODEL_OUTPUT/lat_mig_multcomp_out.RData")
```

Reformat lat migratnt object

```{r}
lmi

lmi$Trait <- "Latitudinal migration\n(Perm. resident vs.\nLat. migrant)"
lmi$Group1 <- c("Perm. resident","Lat. migrant")
lmi$Group <- c("Mature forest","Early successional")


```



Stack all of the data
```{r}
df.effect.sizes <- rbind(lmi,hp,ds,cp,em,c1,c2,fg)
df.effect.sizes$group.x <- with(df.effect.sizes,
                                paste(Trait,Group))
```



Define the order they are to be plotted
```{r}
df.effect.sizes$order <- c(rep("a",2),
                           rep("b",2),
                           rep("c",2),
                           rep("d",2),
                           rep("e",2),
                           rep("f",2),
                           rep("g",2),
                           rep("h",2))
i <- order(df.effect.sizes$order)

df.effect.sizes$Trait <- factor(df.effect.sizes$Trait)

levs <- levels(df.effect.sizes$Trait)[c(8,7,4,3,5,1,2,6)]

df.effect.sizes$Trait <- factor(df.effect.sizes$Trait,
                                levels = rev(levs))

```



## Plot effect sizes 

```{r}


fnt <- 13
pd <- position_dodge(0.395)
plot.out <- ggplot(data = df.effect.sizes,
       aes(y = exp(Effect.size),
           x = Trait,
           #color = Group,
           group = group.x)) +
  geom_errorbar(aes(ymax = exp(Effect.size+ 1.96*SE),
                    ymin = exp(Effect.size- 1.96*SE)),
                width = 0, size = 1.2,
                position = pd) +
    geom_point(size = 5,position = pd,
             aes(color = Group,
                 shape = Group))  +
  theme_bw() + coord_flip() +
  geom_hline(yintercept = 1) + 
  theme(legend.position=c(0.78578954321, 0.935189975)) + 
  theme(panel.grid.minor=element_blank(), 
       panel.grid.major=element_blank()) +
 
  #font size of axes and axes titles 
  theme(axis.title.x = element_text(face="bold", size=fnt+5),
      axis.text.x  = element_text(vjust=0.95, 
      size=fnt,
      face="bold"),
      axis.title.y = element_text(face="bold", 
              size=fnt+5,
            vjust=0.35),
              axis.text.y  = element_text(vjust=0.5, 
              size=fnt,face="bold")) +
  ylab("Annual trend (slope)") +
  xlab("Trait") +
  
#font size of legened
theme(
  #legend.position="right",
  #legend.position=c(0.15, .85),
  legend.title = element_text(colour="black", 
                              size=fnt+3, 
                              face="bold"),
  legend.text = element_text(colour="black", 
                             size = fnt+3, 
                             face = "bold"))
```


Save plot
```{r}
ggsave(plot = plot.out,
       filename = "Figure3.png",width = 8)
```
















# Regression equations
Extract betas and SEs from a model from the list
```{r}
hab.prime.beta <- get.betas(hab.prime.mod.list)
dist.sense.beta <- get.betas(dist.sense.mod.list) #sense2
cons.prior.beta <- get.betas(cons.priority.mod.list)#cons.priority
elev.mig.beta <- get.betas(elev.mig.mod.list) #mig.elevYN
canopy.ob.YN.beta <- get.betas(canopy.ob.YN.mod.list) #canopy.ob.YN
canopy.use.YN.beta <- get.betas(canopy.use.YN.mod.list)#canopy.use.YN
forag.guilde.2lev.beta <- get.betas(forage.guilde.2lev.mod.list)  #forage.guilde.2lev
forage.guilde.3lev.beta <- get.betas(forage.guilde.3lev.mod.list) #forage.guilde2
hab.breadth.beta <- get.betas(hab.breadth.mod.list) #hab.breadth


```







# Habitat breadth data processing

B/c habitat breadth is treated as a continous covariate (it can take on discrete values from 1 to 7) it doesn't jive with the normal workflow.  I therefore extract the model info by hand to put it into the paper.

## Re-load habitat breadth data

```{r}
load("./MODEL_OUTPUT/habitat_breadt_output_list.RData")

str(hab.breadth.output.list,1)

hab.breadth.output.list$anova.table

str(hab.breadth.output.list$all.models,1)

summary(hab.breadth.output.list$all.models$m7a.best)


```
The results are:

beta Year*Habitat breadth = -0.035, SE = 0.009, p < 0.00001


Fixed effects:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)      0.427935   0.295614   1.448 0.147725    
trait           -0.004703   0.084270  -0.056 0.955496    
monthJan        -0.190077   0.129384  -1.469 0.141809    
year.cent        0.120256   0.032367   3.715 0.000203 ***
trait:year.cent -0.035091   0.008837  -3.971 7.16e-05 ***























# Save output
```{r}

date.str <- format(Sys.time(), "%a%b%d%Y")

workspace.str <- paste("workspace_all_models_run",date.str,
                       ".RData", 
                       sep = "")

save.image(file = workspace.str)
```







